{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project 1: Image Feature Extraction and Matching\n",
    "\n",
    "This is Project 1 for [UW CSE P576 Computer Vision](https://courses.cs.washington.edu/courses/csep576/20sp). \n",
    "\n",
    "**Getting Started:** To get started, **[download the source files here (Projects 1 and 2)](https://courses.cs.washington.edu/courses/csep576/20sp/projects/project12/project12.zip \"Project 1 and 2 Source Files\")**. To run the project locally you will need IPython/Jupyter installed, e.g., see instructions at http://jupyter.org/install.html. The notebooks are written for Python 3.x. Launch Jupyter and open `Project1.ipynb`. Alternatively, you can import the standalone version of the notebook into [Colaboratory](https://colab.research.google.com \"Colab\") and run it without installing anything. Use File->Upload Notebook in Colab and open the notebook in `standalone/Project1s.ipynb`.\n",
    "\n",
    "**This project:** In this project you will build an image feature matcher, starting with simple convolution operations and working through interest point detection and descriptor extraction. Once you have a basic feature matcher working, try out some improvements and document your results. If you’re not already familiar with python/numpy, it is recommended to do an introduction such as: http://cs231n.github.io/python-numpy-tutorial. \n",
    "\n",
    "**What to turn in:** Turn in your completed ipynb notebook as well as any source .py files that you modified. Clearly describe any enhancements or experiments you tried in your ipynb notebook. Put everything in a single flat zipfile and upload via the link in Canvas.\n",
    "\n",
    "`version 040820`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os.path\n",
    "from time import time\n",
    "import types\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#import im_util\n",
    "#import interest_point\n",
    "\n",
    "%matplotlib inline\n",
    "# edit this line to change the figure size\n",
    "plt.rcParams['figure.figsize'] = (16.0, 10.0)\n",
    "# force auto-reload of import modules before running code \n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget -nc https://courses.cs.washington.edu/courses/csep576/18sp/projects/project12/pano_images.zip && unzip -n -d data pano_images.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### im_util.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright 2017 Google Inc.\n",
    "\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License.\n",
    "\n",
    "import numpy as np\n",
    "import PIL.Image as pil\n",
    "import scipy.signal as sps\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.ndimage import map_coordinates\n",
    "\n",
    "def convolve_1d(x, k):\n",
    "  \"\"\"\n",
    "  Convolve vector x with kernel k\n",
    "\n",
    "  Inputs: x=input vector (Nx)\n",
    "          k=input kernel (Nk)\n",
    "\n",
    "  Outputs: y=output vector (Nx)\n",
    "  \"\"\"\n",
    "  y=np.zeros_like(x)\n",
    "\n",
    "  \"\"\"\n",
    "  *******************************************\n",
    "  *** TODO: write code to perform convolution\n",
    "  *******************************************\n",
    "\n",
    "  The output should be the same size as the input\n",
    "  You can assume zero padding, and an odd-sized kernel\n",
    "  \"\"\"\n",
    "\n",
    "\n",
    "  \"\"\"\n",
    "  *******************************************\n",
    "  \"\"\"\n",
    "\n",
    "  return y\n",
    "\n",
    "def convolve_rows(im, k):\n",
    "  \"\"\"\n",
    "  Convolve image im with kernel k\n",
    "\n",
    "  Inputs: im=input image (H, W, B)\n",
    "          k=1D convolution kernel (N)\n",
    "\n",
    "  Outputs: im_out=output image (H, W, B)\n",
    "  \"\"\"\n",
    "  im_out = np.zeros_like(im)\n",
    "\n",
    "  \"\"\"\n",
    "  *****************************************\n",
    "  *** TODO: write code to convolve an image\n",
    "  *****************************************\n",
    "\n",
    "  Convolve the rows of image im with kernel k\n",
    "  The output should be the same size as the input\n",
    "  You can assume zero padding, and an odd-sized kernel\n",
    "  \"\"\"\n",
    "\n",
    "\n",
    "  \"\"\"\n",
    "  *****************************************\n",
    "  \"\"\"\n",
    "\n",
    "  return im_out\n",
    "\n",
    "def gauss_kernel(sigma):\n",
    "  \"\"\"\n",
    "  1D Gauss kernel of standard deviation sigma\n",
    "  \"\"\"\n",
    "  l = int(np.ceil(2 * sigma))\n",
    "  x = np.linspace(-l, l, 2*l+1)\n",
    "\n",
    "  # FORNOW\n",
    "  gx = np.zeros_like(x)\n",
    "\n",
    "  \"\"\"\n",
    "  *******************************************\n",
    "  *** TODO: compute gaussian kernel at each x\n",
    "  *******************************************\n",
    "  \"\"\"\n",
    "\n",
    "\n",
    "  \"\"\"\n",
    "  *******************************************\n",
    "  \"\"\"\n",
    "\n",
    "  gx = np.expand_dims(gx,0)\n",
    "  return gx\n",
    "\n",
    "def convolve_gaussian(im, sigma):\n",
    "  \"\"\"\n",
    "  2D gaussian convolution\n",
    "  \"\"\"\n",
    "  imc=np.zeros_like(im)\n",
    "\n",
    "  \"\"\"\n",
    "  ***************************************\n",
    "  *** TODO separable gaussian convolution\n",
    "  ***************************************\n",
    "  \"\"\"\n",
    "\n",
    "\n",
    "  \"\"\"\n",
    "  ***************************************\n",
    "  \"\"\"\n",
    "  return imc\n",
    "\n",
    "def compute_gradients(img):\n",
    "\n",
    "  Ix=np.zeros_like(img)\n",
    "  Iy=np.zeros_like(img)\n",
    "\n",
    "  \"\"\"\n",
    "  ***********************************************\n",
    "  *** TODO: write code to compute image gradients\n",
    "  ***********************************************\n",
    "  \"\"\"\n",
    "\n",
    "\n",
    "  \"\"\"\n",
    "  ***********************************************\n",
    "  \"\"\"\n",
    "  return Ix, Iy\n",
    "\n",
    "def image_open(filename):\n",
    "  \"\"\"\n",
    "  Returns a numpy float image with values in the range (0,1)\n",
    "  \"\"\"\n",
    "  pil_im = pil.open(filename)\n",
    "  im_np = np.array(pil_im).astype(np.float32)\n",
    "  im_np /= 255.0\n",
    "  return im_np\n",
    "\n",
    "def image_save(im_np, filename):\n",
    "  \"\"\"\n",
    "  Saves a numpy float image to file\n",
    "  \"\"\"\n",
    "  if (len(im_np.shape)==2):\n",
    "    im_np = np.expand_dims(im_np, 2)\n",
    "  if (im_np.shape[2]==1):\n",
    "    im_np= np.repeat(im_np, 3, axis=2)\n",
    "  im_np = np.maximum(0.0, np.minimum(im_np, 1.0))\n",
    "  pil_im = pil.fromarray((im_np*255).astype(np.uint8))\n",
    "  pil_im.save(filename)\n",
    "\n",
    "def image_figure(im, dpi=100):\n",
    "  \"\"\"\n",
    "  Creates a matplotlib figure around an image,\n",
    "  useful for writing to file with savefig()\n",
    "  \"\"\"\n",
    "  H,W,_=im.shape\n",
    "  fig=plt.figure()\n",
    "  fig.set_size_inches(W/dpi, H/dpi)\n",
    "  ax=fig.add_axes([0,0,1,1])\n",
    "  ax.imshow(im)\n",
    "  return fig, ax\n",
    "\n",
    "def plot_two_images(im1, im2):\n",
    "  \"\"\"\n",
    "  Plot two images and return axis handles\n",
    "  \"\"\"\n",
    "  ax1=plt.subplot(1,2,1)\n",
    "  plt.imshow(im1)\n",
    "  plt.axis('off')\n",
    "  ax2=plt.subplot(1,2,2)\n",
    "  plt.imshow(im2)\n",
    "  plt.axis('off')\n",
    "  return ax1, ax2\n",
    "\n",
    "def normalise_01(im):\n",
    "  \"\"\"\n",
    "  Normalise image to the range (0,1)\n",
    "  \"\"\"\n",
    "  mx = im.max()\n",
    "  mn = im.min()\n",
    "  den = mx-mn\n",
    "  small_val = 1e-9\n",
    "  if (den < small_val):\n",
    "    print('image normalise_01 -- divisor is very small')\n",
    "    den = small_val\n",
    "  return (im-mn)/den\n",
    "\n",
    "def grey_to_rgb(img):\n",
    "  \"\"\"\n",
    "  Convert greyscale to rgb image\n",
    "  \"\"\"\n",
    "  if (len(img.shape)==2):\n",
    "    img = np.expand_dims(img, 2)\n",
    "\n",
    "  img3 = np.repeat(img, 3, 2)\n",
    "  return img3\n",
    "\n",
    "def disc_mask(l):\n",
    "  \"\"\"\n",
    "  Create a binary cirular mask of radius l\n",
    "  \"\"\"\n",
    "  sz = 2 * l + 1\n",
    "  m = np.zeros((sz,sz))\n",
    "  x = np.linspace(-l,l,2*l+1)/l\n",
    "  x = np.expand_dims(x, 1)\n",
    "  m = x**2\n",
    "  m = m + m.T\n",
    "  m = m<1\n",
    "  m = np.expand_dims(m, 2)\n",
    "  return m\n",
    "\n",
    "def convolve(im, kernel):\n",
    "  \"\"\"\n",
    "  Wrapper for scipy convolution function\n",
    "  This implements a general 2D convolution of image im with kernel\n",
    "  Note that strictly speaking this is correlation not convolution\n",
    "\n",
    "  Inputs: im=input image (H, W, B) or (H, W)\n",
    "          kernel=kernel (kH, kW)\n",
    "\n",
    "  Outputs: imc=output image (H, W, B)\n",
    "  \"\"\"\n",
    "  if (len(im.shape)==2):\n",
    "    im = np.expand_dims(im, 2)\n",
    "  H, W, B = im.shape\n",
    "  imc = np.zeros((H, W, B))\n",
    "  for band in range(B):\n",
    "    imc[:, :, band] = sps.correlate2d(im[:, :, band], kernel, mode='same')\n",
    "  return imc\n",
    "\n",
    "def coordinate_image(num_rows,num_cols,r0,r1,c0,c1):\n",
    "  \"\"\"\n",
    "  Creates an image size num_rows, num_cols\n",
    "  with coordinates linearly spaced in from r0->r1 and c0->c1\n",
    "  \"\"\"\n",
    "  rval=np.linspace(r0,r1,num_rows)\n",
    "  cval=np.linspace(c0,c1,num_cols)\n",
    "  c,r=np.meshgrid(cval,rval)\n",
    "  M = np.stack([r,c,np.ones(r.shape)],-1)\n",
    "  return M\n",
    "\n",
    "def transform_coordinates(coord_image, M):\n",
    "  \"\"\"\n",
    "  Transform an image containing row,col,1 coordinates by matrix M\n",
    "  \"\"\"\n",
    "  M=np.expand_dims(M,2)\n",
    "  uh=np.dot(coord_image,M.T)\n",
    "  uh=uh[:, :, 0, :]\n",
    "  uh=uh/np.expand_dims(uh[:, :, 2],2)\n",
    "  return uh\n",
    "\n",
    "def warp_image(im, coords):\n",
    "  \"\"\"\n",
    "  Warp image im using row,col,1 image coords\n",
    "  \"\"\"\n",
    "  im_rows,im_cols,im_bands=im.shape\n",
    "  warp_rows,warp_cols,_=coords.shape\n",
    "  map_coords=np.zeros((3,warp_rows,warp_cols,im_bands))\n",
    "  for b in range(im_bands):\n",
    "    map_coords[0,:,:,b]=coords[:,:,0]\n",
    "    map_coords[1,:,:,b]=coords[:,:,1]\n",
    "    map_coords[2,:,:,b]=b\n",
    "  warp_im = map_coordinates(im, map_coords, order=1)\n",
    "  return warp_im\n",
    "\n",
    "# allow accessing these functions by im_util.*\n",
    "im_util=types.SimpleNamespace()\n",
    "im_util.convolve_1d=convolve_1d\n",
    "im_util.convolve_rows=convolve_rows\n",
    "im_util.gauss_kernel=gauss_kernel\n",
    "im_util.convolve_gaussian=convolve_gaussian\n",
    "im_util.compute_gradients=compute_gradients\n",
    "im_util.image_open=image_open\n",
    "im_util.image_save=image_save\n",
    "im_util.image_figure=image_figure\n",
    "im_util.plot_two_images=plot_two_images\n",
    "im_util.normalise_01=normalise_01\n",
    "im_util.grey_to_rgb=grey_to_rgb\n",
    "im_util.disc_mask=disc_mask\n",
    "im_util.convolve=convolve\n",
    "im_util.coordinate_image=coordinate_image\n",
    "im_util.transform_coordinates=transform_coordinates\n",
    "im_util.warp_image=warp_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### interest_point.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright 2017 Google Inc.\n",
    "\n",
    "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License.\n",
    "\n",
    "import numpy as np\n",
    "import scipy.ndimage.filters as filters\n",
    "from scipy.ndimage import map_coordinates\n",
    "from matplotlib.patches import Circle\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#import im_util\n",
    "\n",
    "class InterestPointExtractor:\n",
    "  \"\"\"\n",
    "  Class to extract interest points from an image\n",
    "  \"\"\"\n",
    "  def __init__(self):\n",
    "    self.params={}\n",
    "    self.params['border_pixels']=10\n",
    "    self.params['strength_threshold_percentile']=95\n",
    "    self.params['supression_radius_frac']=0.01\n",
    "\n",
    "  def find_interest_points(self, img):\n",
    "    \"\"\"\n",
    "    Find interest points in greyscale image img\n",
    "\n",
    "    Inputs: img=greyscale input image (H, W, 1)\n",
    "\n",
    "    Outputs: ip=interest points of shape (2, N)\n",
    "    \"\"\"\n",
    "    ip_fun = self.corner_function(img)\n",
    "    row, col = self.find_local_maxima(ip_fun)\n",
    "\n",
    "    ip = np.stack((row,col))\n",
    "    return ip\n",
    "\n",
    "  def corner_function(self, img):\n",
    "    \"\"\"\n",
    "    Compute corner strength function in image im\n",
    "\n",
    "    Inputs: img=grayscale input image (H, W, 1)\n",
    "\n",
    "    Outputs: ip_fun=interest point strength function (H, W, 1)\n",
    "    \"\"\"\n",
    "\n",
    "    H, W, _ = img.shape\n",
    "\n",
    "    # FORNOW: random interest point function\n",
    "    ip_fun = np.random.randn(H, W, 1)\n",
    "\n",
    "    \"\"\"\n",
    "    **********************************************************\n",
    "    *** TODO: write code to compute a corner strength function\n",
    "    **********************************************************\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    \"\"\"\n",
    "    **********************************************************\n",
    "    \"\"\"\n",
    "\n",
    "    return ip_fun\n",
    "\n",
    "  def find_local_maxima(self, ip_fun):\n",
    "    \"\"\"\n",
    "    Find local maxima in interest point strength function\n",
    "\n",
    "    Inputs: ip_fun=corner strength function (H, W, 1)\n",
    "\n",
    "    Outputs: row,col=coordinates of interest points\n",
    "    \"\"\"\n",
    "\n",
    "    H, W, _ = ip_fun.shape\n",
    "\n",
    "    # radius for non-maximal suppression\n",
    "    suppression_radius_pixels = int(self.params['supression_radius_frac']*max(H, W))\n",
    "\n",
    "    # minimum of strength function for corners\n",
    "    strength_threshold=np.percentile(ip_fun, self.params['strength_threshold_percentile'])\n",
    "\n",
    "    # don't return interest points within border_pixels of edge\n",
    "    border_pixels = self.params['border_pixels']\n",
    "\n",
    "    # row and column coordinates of interest points\n",
    "    row = []\n",
    "    col = []\n",
    "\n",
    "    # FORNOW: random row and column coordinates\n",
    "    row = np.random.randint(0,H,100)\n",
    "    col = np.random.randint(0,W,100)\n",
    "\n",
    "    \"\"\"\n",
    "    ***************************************************\n",
    "    *** TODO: write code to find local maxima in ip_fun\n",
    "    ***************************************************\n",
    "\n",
    "    Hint: try scipy filters.maximum_filter with im_util.disc_mask\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    \"\"\"\n",
    "    ***************************************************\n",
    "    \"\"\"\n",
    "\n",
    "    return row, col\n",
    "\n",
    "class DescriptorExtractor:\n",
    "  \"\"\"\n",
    "  Extract descriptors around interest points\n",
    "  \"\"\"\n",
    "  def __init__(self):\n",
    "    self.params={}\n",
    "    self.params['patch_size']=8\n",
    "    self.params['ratio_threshold']=1.0\n",
    "\n",
    "  def get_descriptors(self, img, ip):\n",
    "    \"\"\"\n",
    "    Extact descriptors from grayscale image img at interest points ip\n",
    "\n",
    "    Inputs: img=grayscale input image (H, W, 1)\n",
    "            ip=interest point coordinates (2, N)\n",
    "\n",
    "    Returns: descriptors=vectorized descriptors (N, num_dims)\n",
    "    \"\"\"\n",
    "    patch_size=self.params['patch_size']\n",
    "    patch_size_div2=int(patch_size/2)\n",
    "    num_dims=patch_size**2\n",
    "\n",
    "    H,W,_=img.shape\n",
    "    num_ip=ip.shape[1]\n",
    "    descriptors=np.zeros((num_ip,num_dims))\n",
    "\n",
    "\n",
    "    for i in range(num_ip):\n",
    "      row=ip[0,i]\n",
    "      col=ip[1,i]\n",
    "\n",
    "      # FORNOW: random image patch\n",
    "      patch=np.random.randn(patch_size,patch_size)\n",
    "\n",
    "      \"\"\"\n",
    "      ******************************************************\n",
    "      *** TODO: write code to extract descriptor at row, col\n",
    "      ******************************************************\n",
    "      \"\"\"\n",
    "\n",
    "\n",
    "      \"\"\"\n",
    "      ******************************************************\n",
    "      \"\"\"\n",
    "\n",
    "      descriptors[i, :]=np.reshape(patch,num_dims)\n",
    "\n",
    "    # normalise descriptors to 0 mean, unit length\n",
    "    mn=np.mean(descriptors,1,keepdims=True)\n",
    "    sd=np.std(descriptors,1,keepdims=True)\n",
    "    small_val = 1e-6\n",
    "    descriptors = (descriptors-mn)/(sd+small_val)\n",
    "\n",
    "    return descriptors\n",
    "\n",
    "  def compute_distances(self, desc1, desc2):\n",
    "    \"\"\"\n",
    "    Compute distances between descriptors\n",
    "\n",
    "    Inputs: desc1=descriptor array (N1, num_dims)\n",
    "            desc2=descriptor array (N2, num_dims)\n",
    "\n",
    "    Returns: dists=array of distances (N1,N2)\n",
    "    \"\"\"\n",
    "    N1,num_dims=desc1.shape\n",
    "    N2,num_dims=desc2.shape\n",
    "\n",
    "    ATB=np.dot(desc1,desc2.T)\n",
    "    AA=np.sum(desc1*desc1,1)\n",
    "    BB=np.sum(desc2*desc2,1)\n",
    "\n",
    "    dists=-2*ATB+np.expand_dims(AA,1)+BB\n",
    "\n",
    "    return dists\n",
    "\n",
    "  def match_descriptors(self, desc1, desc2):\n",
    "    \"\"\"\n",
    "    Find nearest neighbour matches between descriptors\n",
    "\n",
    "    Inputs: desc1=descriptor array (N1, num_dims)\n",
    "            desc2=descriptor array (N2, num_dims)\n",
    "\n",
    "    Returns: match_idx=nearest neighbour index for each desc1 (N1)\n",
    "    \"\"\"\n",
    "    dists=self.compute_distances(desc1, desc2)\n",
    "\n",
    "    match_idx=np.argmin(dists,1)\n",
    "\n",
    "    return match_idx\n",
    "\n",
    "  def match_ratio_test(self, desc1, desc2):\n",
    "    \"\"\"\n",
    "    Find nearest neighbour matches between descriptors\n",
    "    and perform ratio test\n",
    "\n",
    "    Inputs: desc1=descriptor array (N1, num_dims)\n",
    "            desc2=descriptor array (N2, num_dims)\n",
    "\n",
    "    Returns: match_idx=nearest neighbour inded for each desc1 (N1)\n",
    "             ratio_pass=whether each match passes ratio test (N1)\n",
    "    \"\"\"\n",
    "    N1,num_dims=desc1.shape\n",
    "\n",
    "    dists=self.compute_distances(desc1, desc2)\n",
    "\n",
    "    sort_idx=np.argsort(dists,1)\n",
    "\n",
    "    #match_idx=np.argmin(dists,1)\n",
    "    match_idx=sort_idx[:,0]\n",
    "\n",
    "    d1NN=dists[np.arange(0,N1),sort_idx[:,0]]\n",
    "    d2NN=dists[np.arange(0,N1),sort_idx[:,1]]\n",
    "\n",
    "    ratio_threshold=self.params['ratio_threshold']\n",
    "    ratio_pass=(d1NN<ratio_threshold*d2NN)\n",
    "\n",
    "    return match_idx,ratio_pass\n",
    "\n",
    "def draw_interest_points_ax(ip, ax):\n",
    "  \"\"\"\n",
    "  Draw interest points ip on axis ax\n",
    "  \"\"\"\n",
    "  for row,col in zip(ip[0,:],ip[1,:]):\n",
    "    circ1 = Circle((col,row), 5)\n",
    "    circ1.set_color('black')\n",
    "    circ2 = Circle((col,row), 3)\n",
    "    circ2.set_color('white')\n",
    "    ax.add_patch(circ1)\n",
    "    ax.add_patch(circ2)\n",
    "\n",
    "def draw_interest_points_file(im, ip, filename):\n",
    "  \"\"\"\n",
    "  Draw interest points ip on image im and save to filename\n",
    "  \"\"\"\n",
    "  fig,ax = im_util.image_figure(im)\n",
    "  draw_interest_points_ax(ip, ax)\n",
    "  fig.savefig(filename)\n",
    "  plt.close(fig)\n",
    "\n",
    "def draw_matches_ax(ip1, ipm, ax1, ax2):\n",
    "  \"\"\"\n",
    "  Draw matches ip1, ipm on axes ax1, ax2\n",
    "  \"\"\"\n",
    "  for r1,c1,r2,c2 in zip(ip1[0,:], ip1[1,:], ipm[0,:], ipm[1,:]):\n",
    "    rand_colour=np.random.rand(3,)\n",
    "\n",
    "    circ1 = Circle((c1,r1), 5)\n",
    "    circ1.set_color('black')\n",
    "    circ2 = Circle((c1,r1), 3)\n",
    "    circ2.set_color(rand_colour)\n",
    "    ax1.add_patch(circ1)\n",
    "    ax1.add_patch(circ2)\n",
    "\n",
    "    circ3 = Circle((c2,r2), 5)\n",
    "    circ3.set_color('black')\n",
    "    circ4 = Circle((c2,r2), 3)\n",
    "    circ4.set_color(rand_colour)\n",
    "    ax2.add_patch(circ3)\n",
    "    ax2.add_patch(circ4)\n",
    "\n",
    "def draw_matches_file(im1, im2, ip1, ipm, filename):\n",
    "  \"\"\"\n",
    "  Draw matches ip1, ipm on images im1, im2 and save to filename\n",
    "  \"\"\"\n",
    "  H1,W1,B1=im1.shape\n",
    "  H2,W2,B2=im2.shape\n",
    "\n",
    "  im3 = np.zeros((max(H1,H2),W1+W2,3))\n",
    "  im3[0:H1,0:W1,:]=im1\n",
    "  im3[0:H2,W1:(W1+W2),:]=im2\n",
    "\n",
    "  fig,ax = im_util.image_figure(im3)\n",
    "  col_offset=W1\n",
    "\n",
    "  for r1,c1,r2,c2 in zip(ip1[0,:], ip1[1,:], ipm[0,:], ipm[1,:]):\n",
    "    rand_colour=np.random.rand(3,)\n",
    "\n",
    "    circ1 = Circle((c1,r1), 5)\n",
    "    circ1.set_color('black')\n",
    "    circ2 = Circle((c1,r1), 3)\n",
    "    circ2.set_color(rand_colour)\n",
    "    ax.add_patch(circ1)\n",
    "    ax.add_patch(circ2)\n",
    "\n",
    "    circ3 = Circle((c2+col_offset,r2), 5)\n",
    "    circ3.set_color('black')\n",
    "    circ4 = Circle((c2+col_offset,r2), 3)\n",
    "    circ4.set_color(rand_colour)\n",
    "    ax.add_patch(circ3)\n",
    "    ax.add_patch(circ4)\n",
    "\n",
    "  fig.savefig(filename)\n",
    "  plt.close(fig)\n",
    "\n",
    "def plot_descriptors(desc,plt):\n",
    "  \"\"\"\n",
    "  Plot a random set of descriptor patches\n",
    "  \"\"\"\n",
    "  num_ip,num_dims = desc.shape\n",
    "  patch_size = int(np.sqrt(num_dims))\n",
    "\n",
    "  N1,N2=2,8\n",
    "  figsize0=plt.rcParams['figure.figsize']\n",
    "  plt.rcParams['figure.figsize'] = (16.0, 4.0)\n",
    "  for i in range(N1):\n",
    "    for j in range(N2):\n",
    "      ax=plt.subplot(N1,N2,i*N2+j+1)\n",
    "      rnd=np.random.randint(0,num_ip)\n",
    "      desc_im=np.reshape(desc[rnd,:],(patch_size,patch_size))\n",
    "      plt.imshow(im_util.grey_to_rgb(im_util.normalise_01(desc_im)))\n",
    "      plt.axis('off')\n",
    "\n",
    "  plt.rcParams['figure.figsize']=figsize0\n",
    "\n",
    "def plot_matching_descriptors(desc1,desc2,desc1_id,desc2_id,plt):\n",
    "  \"\"\"\n",
    "  Plot a random set of matching descriptor patches\n",
    "  \"\"\"\n",
    "  num_inliers=desc1_id.size\n",
    "  num_ip,num_dims = desc1.shape\n",
    "  patch_size=int(np.sqrt(num_dims))\n",
    "\n",
    "  figsize0=plt.rcParams['figure.figsize']\n",
    "\n",
    "  N1,N2=1,8\n",
    "  plt.rcParams['figure.figsize'] = (16.0, N1*4.0)\n",
    "\n",
    "  for i in range(N1):\n",
    "    for j in range(N2):\n",
    "      rnd=np.random.randint(0,num_inliers)\n",
    "\n",
    "      desc1_rnd=desc1_id[rnd]\n",
    "      desc2_rnd=desc2_id[rnd]\n",
    "\n",
    "      desc1_im=np.reshape(desc1[desc1_rnd,:],(patch_size,patch_size))\n",
    "      desc2_im=np.reshape(desc2[desc2_rnd,:],(patch_size,patch_size))\n",
    "\n",
    "      ax=plt.subplot(2*N1,N2,2*i*N2+j+1)\n",
    "      plt.imshow(im_util.grey_to_rgb(im_util.normalise_01(desc1_im)))\n",
    "      plt.axis('off')\n",
    "      ax=plt.subplot(2*N1,N2,2*i*N2+N2+j+1)\n",
    "      plt.imshow(im_util.grey_to_rgb(im_util.normalise_01(desc2_im)))\n",
    "      plt.axis('off')\n",
    "\n",
    "  plt.rcParams['figure.figsize'] = figsize0\n",
    "\n",
    "# allow accessing these functions by interest_point.*\n",
    "interest_point=types.SimpleNamespace()\n",
    "interest_point.InterestPointExtractor=InterestPointExtractor\n",
    "interest_point.DescriptorExtractor=DescriptorExtractor\n",
    "interest_point.draw_interest_points_ax=draw_interest_points_ax\n",
    "interest_point.draw_interest_points_file=draw_interest_points_file\n",
    "interest_point.draw_matches_ax=draw_matches_ax\n",
    "interest_point.draw_matches_file=draw_matches_file\n",
    "interest_point.plot_descriptors=plot_descriptors\n",
    "interest_point.plot_matching_descriptors=plot_matching_descriptors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution and Image Filtering [25%]\n",
    "\n",
    "Start by writing code to perform convolution in 1D. Open `im_util.py` and edit the function `convolve_1d`. You should use only basic numpy array operations and loops. Don't worry about efficiency for now. You should see small errors compared to the reference numpy version.\n",
    "\n",
    "Note that convolution and correlation are the same under a simple manipulation of the kernel (what is it?). For what kernels are convolution and correlation results identical?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Test of convolve_1d\n",
    "\"\"\"\n",
    "print('[ Test convolve_1d ]')\n",
    "x = (np.random.rand(20)>0.8).astype(np.float32)\n",
    "k = np.array([1,3,1])\n",
    "y1 = im_util.convolve_1d(x, k)\n",
    "y2 = np.convolve(x, k, 'same')\n",
    "y3 = np.correlate(x, k, 'same')\n",
    "print(' convolve error = ', np.sum((y1-y2)**2))\n",
    "print(' correlate error = ', np.sum((y1-y3)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now convolve a 2D image with a 1D kernel. Before you begin, get some image data by running `get_data.sh` in the `data` directory. Then complete the function `convolve_rows` in `im_util.py` by convolving every row of the image by the kernel. Run the code below and check that the image output is sensible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Test of convolve_image\n",
    "\"\"\"\n",
    "image_filename='data/test/100-0038_img.jpg'\n",
    "\n",
    "print('[ Test convolve_image ]')\n",
    "im = im_util.image_open(image_filename)\n",
    "k = np.array([1,2,3,4,5,6,5,4,3,2,1])\n",
    "print(' convolve_rows')\n",
    "t0=time()\n",
    "im1 = im_util.convolve_rows(im, k)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "print(' scipy convolve')\n",
    "t0=time()\n",
    "im2 = im_util.convolve(im, np.expand_dims(k,0))\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "print(' convolve_image error =', np.sum((im1-im2)**2))\n",
    "\n",
    "# optionally plot images for debugging\n",
    "#im1_norm=im_util.normalise_01(im1)\n",
    "#im2_norm=im_util.normalise_01(im2)\n",
    "#ax1,ax2=im_util.plot_two_images(im1_norm, im2_norm)          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will probably find that the scipy convolve runs much faster than your version. To speed things up you can use this version (`im_util.convolve`) for all subsequent experiments. Note that this performs a general 2D convolution with a 2D kernel as input. \n",
    "\n",
    "Now write code to perform Gaussian blurring. First implement the function `gauss_kernel` to compute a 1D Gaussian kernel. Then complete `convolve_gaussian` to perform a separable convolution with this kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Gaussian blurring test\n",
    "\"\"\"\n",
    "print('[ Test convolve_gaussian ]')\n",
    "\n",
    "sigma=4.0\n",
    "k=im_util.gauss_kernel(sigma)\n",
    "print(' gauss kernel = ')\n",
    "print(k)\n",
    "\n",
    "t0=time()\n",
    "im1 = im_util.convolve_gaussian(im, sigma)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "ax1,ax2=im_util.plot_two_images(im, im1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now write code to compute horizontal and vertical gradients in the function `compute_gradients`. Use an explicit kernel that is convolved in each direction (i.e., do not use a built-in function such as `numpy.gradient`). Run the code below and check that the outputs look sensible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Gradient computation test\n",
    "\"\"\"\n",
    "print('[ Test gradient computation ]')\n",
    "img = np.mean(im,2,keepdims=True)\n",
    "Ix,Iy = im_util.compute_gradients(img)\n",
    "\n",
    "# copy greyvalue to RGB channels\n",
    "Ix_out = im_util.grey_to_rgb(im_util.normalise_01(Ix))\n",
    "Iy_out = im_util.grey_to_rgb(im_util.normalise_01(Iy))\n",
    "\n",
    "im_util.plot_two_images(Ix_out, Iy_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interest Point Extractor [25%]\n",
    "\n",
    "Now you will use these convolution functions to implement a corner or interest point detector. Choose a well known detector, such as Harris or DoG, and implement the interest point strength function in `corner_function` of `interest_point.py`. Run the code below to visualise your corner function output. Next detect corners as local maxima of this function by filling in `find_local_maxima` in the same file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compute corner strength function\n",
    "\"\"\"\n",
    "print('[ Compute corner strength ]')\n",
    "ip_ex = interest_point.InterestPointExtractor()\n",
    "ip_fun = ip_ex.corner_function(img)\n",
    "\n",
    "# normalise for display\n",
    "[mn,mx]=np.percentile(ip_fun,[5,95])\n",
    "small_val=1e-9\n",
    "ip_fun_norm=(ip_fun-mn)/(mx-mn+small_val)\n",
    "ip_fun_norm=np.maximum(np.minimum(ip_fun_norm,1.0),0.0)\n",
    "\n",
    "\"\"\"\n",
    "Find local maxima of corner strength\n",
    "\"\"\"\n",
    "print('[ Find local maxima ]')\n",
    "row, col = ip_ex.find_local_maxima(ip_fun)\n",
    "ip = np.stack((row,col))\n",
    "\n",
    "ax1,ax2=im_util.plot_two_images(im_util.grey_to_rgb(ip_fun_norm),im)\n",
    "interest_point.draw_interest_points_ax(ip, ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Descriptors and Matching [25%]\n",
    "\n",
    "Now let's match our interest points. Start by extracting a very simple descriptor that is simply a patch of pixels around the interest point. To do this, fill in the function `get_descriptors` in `interest_point.py`. The following code outputs a random set of normalised descriptor patches. Check that the output looks sensible. Once you have this working, try varying the sample spacing in your descriptor patch. \"Sample spacing\" here means distance in pixels in the base image between the pixels in your descriptor patch. What problem exists with sample spacings > 1 pixel? How can this be fixed?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Extract descriptors\n",
    "\"\"\"\n",
    "print('[ Extract descriptors ]')\n",
    "desc_ex=interest_point.DescriptorExtractor()\n",
    "descriptors=desc_ex.get_descriptors(img, ip)\n",
    "interest_point.plot_descriptors(descriptors,plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now match descriptors between a pair of images. Run the following two code blocks to extract your interest points and extract and match descriptors. The second code block calls a function to perform nearest-neighbour matching of descriptors and filtering using a ratio test. Take a look at the code and check you understand how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Read a pair of input images and extract interest points\n",
    "\"\"\"\n",
    "image_dir='data/test'\n",
    "#im_filename1=image_dir+'/100-0023_img.jpg'\n",
    "#im_filename2=image_dir+'/100-0024_img.jpg'\n",
    "im_filename1=image_dir+'/100-0038_img.jpg'\n",
    "im_filename2=image_dir+'/100-0039_img.jpg'\n",
    "\n",
    "im1 = im_util.image_open(im_filename1)\n",
    "im2 = im_util.image_open(im_filename2)\n",
    "\n",
    "img1 = np.mean(im1, 2, keepdims=True)\n",
    "img2 = np.mean(im2, 2, keepdims=True)\n",
    "\n",
    "print('[ find interest points ]')\n",
    "t0=time()\n",
    "ip_ex = interest_point.InterestPointExtractor()\n",
    "ip1 = ip_ex.find_interest_points(img1)\n",
    "print(' found '+str(ip1.shape[1])+' in image 1')\n",
    "ip2 = ip_ex.find_interest_points(img2)\n",
    "print(' found '+str(ip2.shape[1])+' in image 2')\n",
    "t1=time()\n",
    "print(' % .2f secs ' % (t1-t0))\n",
    "\n",
    "print('[ drawing interest points ]')\n",
    "ax1,ax2=im_util.plot_two_images(im1,im2)\n",
    "t0=time()\n",
    "interest_point.draw_interest_points_ax(ip1, ax1)\n",
    "interest_point.draw_interest_points_ax(ip2, ax2)\n",
    "t1=time()\n",
    "print(' % .2f secs ' % (t1-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Extract and match descriptors\n",
    "\"\"\"\n",
    "print('[ extract descriptors ]')\n",
    "t0=time()\n",
    "desc_ex = interest_point.DescriptorExtractor()\n",
    "desc1 = desc_ex.get_descriptors(img1, ip1)\n",
    "desc2 = desc_ex.get_descriptors(img2, ip2)\n",
    "t1=time()\n",
    "print(' % .2f secs' % (t1-t0))\n",
    "\n",
    "# Uncomment the following lines to use SIFT descriptors\n",
    "# Note: you'll need to install cyvlfeat, e.g., conda install -c menpo cyvlfeat\n",
    "\n",
    "#from cyvlfeat import sift \n",
    "#frames1,desc1=sift.sift(img1,compute_descriptor=True,n_levels=1)\n",
    "#frames2,desc2=sift.sift(img2,compute_descriptor=True,n_levels=1)\n",
    "#ip1=(frames1.T)[0:2,:]\n",
    "#ip2=(frames2.T)[0:2,:]\n",
    "#desc1=desc1.astype(np.float)\n",
    "#desc2=desc2.astype(np.float)\n",
    "\n",
    "print('[ match descriptors ]')\n",
    "match_idx,ratio_pass=desc_ex.match_ratio_test(desc1, desc2)\n",
    "num_ratio_pass=np.sum(ratio_pass)\n",
    "\n",
    "ipm=ip2[:,match_idx]\n",
    "\n",
    "ip1r=ip1[:,ratio_pass]\n",
    "ip2r=ipm[:,ratio_pass]\n",
    "\n",
    "N1,num_dims=desc1.shape\n",
    "print(' Number of interest points = '+str(N1))\n",
    "print(' Number of matches passing ratio test = '+str(num_ratio_pass))\n",
    "\n",
    "ax1,ax2=im_util.plot_two_images(im1,im2)\n",
    "interest_point.draw_matches_ax(ip1r, ip2r, ax1, ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code visualises matched descriptor patches. Can you distinguish the correct and incorrect matches? (reload to get another random sample)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plot descriptors for matched points\n",
    "\"\"\"\n",
    "interest_point.plot_matching_descriptors(desc1,desc2,np.arange(0,ip1.shape[1]),match_idx,plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing and Improving Feature Matching [25%]\n",
    "\n",
    "Try varying the `ratio_threshold` parameter in the descriptor matcher (`DescriptorExtractor` class params). What are good settings for this parameter? If everything is working, you should see a good set of correctly matched points (aim for about 100 or more). Experiment with your interest point and descriptor implementations to find which other parameters are important and try to get a good set of matches. Try out a new idea of your own to improve interest points or descriptors, and record your findings in the notebook below. Make sure to include sufficient figures/tables and explanations to demonstrate your base case and the effect of your modifications.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### TODO experiments with your detector/descriptors\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
